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Abstract. Composite systems, where couplings are of two types, a combination 
of strong dilute and weak dense couplings of Ising spins, are examined through the 
replica method. The dilute and dense parts are considered to have independent 
canonical disordered or uniform bond distributions; mixing the models by variation of a 
parameter 7 alongside inverse temperature j3 we analyse the respective thermodynamic 
solutions. We describe the variation in high temperature transitions as mixing occurs; 
in the vicinity of these transitions we exactly analyse the competing effects of the dense 
and sparse models. By using the replica symmetric ansatz and population dynamics 
we described the low temperature behaviour of mixed systems. 



1. Introduction 

Understanding phenomena arising in many body systems through mean field analysis of 
simple models has provided insight into many problems in physics, theoretical computer 
science, telecommunication, biology and elsewhere [MPV87, NisOl]. Statistical 
mechanics describes aspects of macroscopic behaviour in interacting systems of many 
elements, and methods originating in the study of spin-glasses (SG) have been extended 
to explore many interesting model disordered systems. These statistical descriptors of 
behaviour often prove to be a sufficient descriptor of all interesting bulk behaviour; in 
some applications, such as channel coding and theoretical computer science, they also 
provide benchmarks even where the large system is relatively small and the randomness 
assumed does not quite match the true system conditions. Moreover, methods developed 
within the statistical physics community gave rise to the development of efficient 
inference algorithms widely used in telecommunication, probabilistic modelling and 
theoretical computer science. 

Many of these models consider range-free interactions of a single topological type, 
most commonly nearest neighbour interactions in finite dimensional, fully connected, 
or sparse random graphs. Even in systems not conforming strongly to the particular 
topology, insight into many phenomena can be obtained and the formulations can 
allow for exact analysis. We propose that certain systems may be well described by 
a combination of two canonical topologies - such a property may be phenomenological 
or could be a deliberately engineered feature of a system. 
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Amongst the best understood topologies are those amenable to mean field theory 
models that are infinite dimensional. The canonical mean field model is that of a 
fully connected graph. This model may be analysed exactly for the cases of uniform 
interactions and certain random ensembles, most famously the celebrated Sherrington 
Kirkpatrick (SK) model [SK75]. Simplification of the analysis in the disordered case is 
often possible through noting the ability to describe many combinations of interactions 
by a Gaussian due to the central limit theorem. The sparse mean field model has 
each spin coupled only to a small number of other spins. This creates a topology 
that is also described as infinite dimensional though the model is perceived as more 
representative of most complex systems due to some notion of a neighbourhood being 
maintained. Analysis is simplified by utilising the asymptotic cycle free property (Bethe 
approximation) of certain random graph ensembles. Models which do not allow use of 
the Bethe approximation or central limit theorem are in general difficult to analyse. 

The motivation for studying this system is that it is amongst the simplest composite 
models involving a combination of dilute and dense couplings. We anticipate the 
understanding of these systems to be important in engineering applications. With 
miniaturisation of technology, for example computer chips, the paradigm by which the 
interaction amongst components can be strictly controlled may be invalid. Modelling 
the effect of non-engineered couplings by independent noise may be invalid in some 
scenarios, as it is possible that these additional interactions form a strongly correlated 
network allowing for example non-trivial phase transitions or metastable features. There 
is also the possibility that systems may be engineered deliberately with a combination 
of dilute strong, and weak dense interactions either to make them more robust or to 
exploit specific properties of the composite system. In multi-user channel coding, for 
instance, this may make the communication process more robust against different types 
of noise, de-synchronisation or malicious attacks. Only recently, a special case of the 
system studied here was suggested as a model for studying the resilience of networks 
against attacks [HM]. 

The paper is organised as follows. In section 2 we introduce the model studied and 
basic definitions, followed by an analysis section 3 that explains briefly the derivation 
based on the replica method. In sections 4 and 5 we investigate the high and low 
temperature solutions, respectively, followed by a discussion of the numerical solutions 
obtained via population dynamics section 6. Finally, we will present our conclusions in 
section 7. 

2. The model 

We present analysis for an ensemble of disordered systems of N Ising spins Sj = ±1. 
The couplings consists of a strong part which is non-zero on only a fraction pN/ 2 of the 
possible links, and a weak part which is present on all links, this defines the Hamiltonian 
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consisting of sparse and dense parts 

H(S) = H S (S) + H D {S) + ]T hiSi (1) 

i 

along with an external field. The dense and sparse parts are taken to be 
H D (S) = -J^bibjJ^o-iO-j, , 

(ij) 

H S (S) = Y,h;hj.\„ 4 a, a,. (2) 

(ij) 

where () indicates the set of ordered (distinct) indices throughout the paper. 
Components of the sparse connectivity matrix A take the value 1 if a link exists between 
the corresponding sites and zero otherwise; the coupling matrix J takes random values 
from a given distribution and b determines an artificial disorder (alignment) in the 
coupling strengths, analogous to the Mattis model [Mat]. For our analysis we can take 
b s = 1 since only the relative alignment is influential in determining system properties. 

2.0.1. Mixing of models In order to investigate the combination of these subsystems 
we propose the introduction of two parameters, an inverse temperature /3, and a mixing 
parameter 7 e (0, 1). The mixing parameter acts so that the couplings in the sparse 
part increase monotonically from zero with 7 to their full value, conversely the couplings 
of the dense part decrease monotonically to zero. There is some flexibility in how this 
might be applied, for example if the couplings decrease/increase linearly with 7 one may 
write the Hamiltonian as 

H(S) = -fHs(S) + (1 - -y)H D (S) + £ hiSi . (3) 

i 

This is the simplest composition method we may use but alternative rescalings of the 
couplings may also be sensible. Unattractive features include that the ratio of variance 
to mean (J /J) in coupling strengths is not maintained as 7 varies. 

The mixing of models by rescaling of the couplings is not a unique way to 
consider the composition of two such systems. A sensible alternative valid in some 
composite systems would be a doping one, keeping the dense part constant and gradually 
introducing additional sparse bonds (variation of A) [HM] . 

2.0.2. Definition of the ensemble We consider different ensembles determined by a set 
of parameters X = j J , J, p, 0, h s , h D , x, &| defined as follows. The model consists of 
independently generated dense and sparse subsystems. Dense Couplings are sampled 
independently for each link from a distribution of mean J bibj/N and variance J/N 
and non-divergent higher moments (the scaling in iV is standard [SK75]). The Mattis 
model- like part b describes some non-trivial orientation of the spins with bi — ±1 
sampled independently at each site according to the mean b. For the sparse part we 
have that the connectivity matrix, A, is a sample from an Erdos-Renyi random graph 
ensemble [ER70] of mean connectivity p with couplings sampled independently for each 
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link from a distribution on the real line, 0, with non-divergent moments, and vanishing 
measure on zero. The external fields are sampled for each site from a distribution of 
mean h s + bih D , and variance x 2 , the parameters hs,fiD and x 2 are conjugate variables 
in the free energy to the order parameters for sparse and dense aligned ferromagnetic 
moments, and the spin glass moment. 

The model contains a wide range of parameters which we believe are sufficient to 
describe the mixing of many canonical Ising spin mean field models. In the case of 
7^0 the model reduces to the SK model [SK75] (up to artificial disorder), whereas at 
7 — > 1 the model is Viana Bray (VB) [VB85]. By tuning parameters one can also find 
the ferromagnetic, antiferromagnetic and Mattis models [Mat] , the relevant orientation 
in mixing (6) proves important in determining system properties. Small perturbations 
of the SK and VB models by random Hamiltonians is a subject much studied, especially 
in the context of temperature variation and stochastic stability [BS07, Par, Tal03]. We 
understand that the set of perturbations represented by variation of an infinitesimal 
variation of 7 from or 1 probably falls into the classes which are incapable of changing 
the structure of thermodynamic states - provided we break any interaction symmetries 
by addition of small external fields, and so we expect no transitions in these limits, 
which is both an observed and intuitive assumption. 



3. Replica method and exact analysis 

3.1. Self- averaged free energy calculation 

The analysis is rather standard and is carried out by the replica method, we note only 
important key steps and definitions in this section. Details of the derivation are found 
in Appendix A. Throughout analysis we consider only the leading order contribution, 
in N, to the free energy density, and in all coefficients we assume the asymptotic value 
in n where possible for brevity. We first write the free energy function making use of 
the replica trick [MPV87] . 

/= Hm-i-log £ e X p{-PH(S 1 ,...,S N )} , (4) 

iV->oo pl\ 



1 d n 
= lim --777 lim — T~J 

n^oo j3N n^o dn ^ 



exp{-[3H(S?,...,S a N )} 



Ca col 



We are interested in the behaviour of a particular sample drawn from the ensemble 
described by the parameters X. We anticipate that a sufficient description of any typical 
instance will be given by the free energy averaged over instances of the disorder (self 
averaging assumption). 

After taking the quenched averages of couplings and some manipulation of the 
equation form we are able to describe the self averaged free energy density by 

(/> = jjm -^lm^/d$d$exp{-JVA(n,J,$,$)} , (5) 



On composite systems 5 

where $, $ = |7r(cr), 7r(<r), q a , q a , q{aia 2 )i 3(ma 2 )} are the set of integration variables 
introduced to allow exact site factorisation. These may be invoked as order parameters 
for various phases. Considering only leading order terms in N one can present A is in 
the general composite case of two types of disordered couplings 

A= -£5a9a- E 9<aia 2 >9<aia 2 > ~ E ^"O^* 7 ) - log / ^ exp { ^} \ 



(aia 2 > 



cr 1 ,(T2 \\ ( a ) I 4>{x) 

(3J 
2 



2 + E (^<«ia2>) 2 
(ai«2) 



Efe) 2 -/3E(^E^K + ^)-/3V(|+ E 9<«i<»: 



* = " E ^ - E 9( aia2 )^V Q2 - tt(<x) . (6) 

a («i«2) 

The scaling with 7 is absorbed in the coupling distribution <p,J Q and J. 

From here it is possible to evaluate the integral by the saddle point method 
anticipating the large N limit. The self-consistent equations obtained by the 
extremisation of the exponent may be evaluated exactly at high temperature and by a 
qualified approximation at lower temperature. The saddle point equations for the order 
parameters are 

q a =^(j2bo- a expx) ; q a = f3Joq a fihP (7) 
\<j I b 

5(aia 2 > = M{Y^a ai a a2 exp X \ ; q {aia2) = -f3 2 Jq {aia2) - (3 2 X 2 (8) 

\ a lb 
7r(«r) =M(expX) b ; 



7r(<ri) = -p|>(<r 2 )^exp {/fc5>?a?}^ 



-1 -^E^ 81 ( 9 ) 

^(x) / "1 



AT = (]Texp*\ 
\ a lb 



from which the conjugate order parameters <E> may be eliminated. This representation 
is convenient for making general ansatze [Mon98] on the order parameters in (39). 

3.2. Alternative formulation 

It is useful to represent the exponent (??) by choosing a standard parameterisation 
for Ti allowing components to be associated with different types of order. One can 
equivalently generate this representation directly by choosing different order parameters 
and expansions in the replica calculation, the representation remains a standard 
one [VB85, WS88]. Working at the saddlepoint we are able to substitute saddlepoint 
definitions for n,q a and q{ aia2 ) to obtain an equation in only ir(cr),q a and q{ ai a 2 )- We 
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then reparameterise 7r, the generating function of the order parameters, by the complete 
expansion 

7T(<7) = 1 + E + E ?<«i«2>* ai * aa + E ^...a,)^ " " " > (10) 

" («i« 2 ) (=3 

The components go as well as q a represent spin-spin correlation order parameters 
between replica [NisOl]. In the case 6=1 the distinction between q a and q a is artificial, 
as emerges from the calculation. We can rewrite the exponent at the saddlepoint as 

A = -log(Eexp{*}\ +^Efe) 2 + / ^ PT2 E (?<« 10a >) 2 

\ff I b a (oia 2 ) 

+ ^E(^) 2 + E^ E (9«i...«.) 2 - 

« 2=3 (ai...ai) 

X= J2(3(KJoq a + h D )+pT iqa + h s )o- a + (3 2 £ ((J + pT 2 )q( ai a 2 ) + X 2 )o- ai cr a2 
Q (0:102) 

+ PE^ E (Zax..^ 1 ---^- (11) 

Z=3 (ai...aj) 

up to addition of constants. The quantities % are to leading order in n 

% = J d<j>{x) tanh^/te) (12) 

with being the coupling distribution. The new order parameters must satisfy the set 
of coupled saddlepoint equations 



q {ai ...a p )= (AT^K 1 ...^}expA'\ , 
\ a Ih 



(13) 



with (7) still applying. Non-zero order solutions of many types may exist but solutions 
are non-trivial except at high temperature and zero exteral fields. In the next section 
we examine solutions for vanishing external fields hg = ho = x — emergent in the 
high temperature regime. 



4. High temperature solutions at zero external field 

In the case of zero external fields there exist phases classified by the set of non-zero order 
parameters. For the low (3 regime the unique stable solution to these equations is the 
paramagnetic one (P), with all order parameters zero (13). It is possible to determine 
a hierarchy of critical temperatures [VB85] for the emergence of different types of non- 
paramagnetic order solution by writing the saddle point equations (13) in cases with 
all but one type of order parameter types zero. A single spin (ferromagnetic, F) order 
emerges subject to a non-zero solution of 

( ^ a \ = ( tan HP J oq a ) + btanh(pT 1 q a ) \ 
\ q a J \ 6tanh(/3J g a ) + tanh(pTig a ) j ' 

A Spin-Glass (SG) order requires a non-zero solution to 

= tanh((/3 2 J + pT 2 )g ( 

) • (15) 
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Due to the convexity of the tanh function the criteria is described equivalently by a 
linear expansion 

A 5 = (3 2 J + P T 2 ] > 1 . (16) 

Higher order solutions, indicated by q( ai ... ap ) emerge subject to the criteria pT p > 1, by 
a similar linearisation to (16). However, no solution of this type can emerge at higher 
temperature than that indicated by equation (16) as the following inequality holds for 
arbitrary coupling distribution [VB85] 

(/3 2 J + pT 2 )>A^, (17) 

defining A^ = pT p for p > 2, the inequality becomes strict in the case J > 0. 

One can apply similar arguments in convexity to aid solution finding in the second 
case (14), we will restrict attention to the cases 6=1 and 6 = 0. In the first case the 
right-hand-side of (14) is identical in the two components so that the order emerges 
uniquely along the direction t a oc q a = q a . The criteria for a non-zero solution is 

'X+ = (3 Jo + pTt] > 1 . (18) 

For the case 6 = the situation is also relatively simple, by the convexity properties of 
the tanh it is sufficient to linearise (15) in q a , q a . The independent processes can then 
yield solutions depending on one of two eigenvalues meeting the criteria 

A± = \ (P J o + P r i ± y/(M ~ P?i) 2 + ^JopTx)] >1, (19) 

written to be inclusive of both cases 6 = and 6=1 (18). Except in the case of a critical 
temperature with pT\ = /3J = 1 the largest eigenvalue \ F = A+ will determine the type 
of the emergent one spin order. We use the generalised order parameter t a oc v ±q a + v 2 q a 
to correspond to the mode A+, v depending on the type of Ferromagnetic order. We 
assume throughout the paper that the discrete symmetry in solutions ±t a is broken by 
a small external field (h s or h D ) aligned with the positive solution. When A^ \ f 
one must consider both q a and q a becoming non-zero simultaneously, the consequences 
of generalising the following sections' analysis to include this case are considered in 
Appendix B.l. 

Equations (16), (17) and (19) indicate that for any mixed system (J or J ^ 0) there 
is a transition at some temperature towards a ferromagnetic or SG phase. The effect of p, 
which indicates a percolation in the sparse couplings only if p > 1, generates no obvious 
criticality in the expressions, except in cases where J < (J = 0) then p must exceed 
1 for the possibility of a ferromagnetic (SG) high temperature transition, regardless 
of bond strength or distribution in the sparse part. In this scenario of negative mean 
coupling in the dense subsystem, A^ (19) is a convex, not monotonically increasing, 
function of (5. This means for example that by continuous variation of several ensemble 
parameters one is able to create a first order Paramagnetic-Ferromagnetic (sparsely 
aligned) transition, with the upper critical temperature changing discontinuously in 7. 
Details may depend on the sparse coupling distribution 0. 
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There exist a number of alternative rigorous methods to attain such a set of 
paramagnetic high temperature transition points without the use of replicas. In 
particular we note the results [Ost06], which allow the high temperature transitions 
to be found by transformation of a disordered Ising spin model to a uniform interaction 
Ising model for many topologies, thereby allows easy identification of transition points. 



4-1. High temperature auxiliary system 

We generate a system to describe the non-Paramagnetic phase in the vicinity of the 
high temperature transition with either or both Ai = — 1 or A2 = X SG — 1 small 
and positive, and with A^ < 1 and A^ < 1 (there is only one potentially non-zero 
ferromagnetic alignment). In this scenario we may describe the leading order behaviour 
to second order in A exactly by an auxiliary saddlepoint equation - which may be 
subject to exact analysis. This approach is motivated by related methods for sparse 
SG [VB85, Mot87]. 

We first introduce some important definitions chosen to be compatible with [VB85] 

= Xf-i; r * = h~ l] n>2 = aw" 1; 

t a = A+(wig Q + v 2 q a ) + Pivth 3 + v 2 h D ); Q( ai a 2 ) = X s q( aiCt2 ) + P 2 X 2 5 (20) 

Q(aia 2 a 3 ) A^ ^ '^(ctio^og) 1 Q '(«ia2«3a4) ^ ^ Q{aia 2 a 3 ai) 

Notation () is to indicate any permutation in the set of indices, the new order parameters 
are identical for any ordering of indices. The notation () is dropped for abbreviation 
in the parameters henceforth. The components of v = (t>i,t>2) describe the different 
possible alignments of the ferromagnetic order at the high temperature boundary: (1,1) 
for b — 1 and (1, 0) [(0, 1)] for b = with pTi < [>]/3 J , respectively. 

At a certain temperature, the minimum j3 which satisfies either (16) or (19), a new 
phase continuously emerge from the paramagnetic solution which is either ferromagnetic 
or SG, respectively. A sufficient description of these solutions in the vicinity of this high 
temperature boundary is given by an expansion in the restricted set of order parameters: 
ferromagnetic t a and SG Q ai a 2 under some ansatz. The rescaling of order parameters 
(20) is to abbreviate the expressions, and we consider the case of negligable external fields 
in examining the transitions. For X x < 1 the saddle point solution for the X th order 
parameter is zero unless there is a non-zero component in the higher order parameters, 
we say the order is induced. For example a non-zero t a induces order in Q aia2 , Q ai a 2 a 3 
when X s , A^ < 1, and Q ai a 2 induces order in Q ai a 2 a 3 a 4 when A^ < 1. We include order 
parameters upto A th order, this set is sufficient to describe the leading order behaviours 
of the phases. These inclusions discriminate the approach from a simple mean field 
(fully connected) approximation to the correlation structure of the sparse subsystem. 

Calculations of the significant higher order parameters through the saddle point 
equations (13) and of Tr(expA') are undertaken in Appendix B. The auxiliary 
expression for A (??) can then be found by expansion in t a ,Q aia2 ,Q aia2a3 ,Q aia2a3a4 
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to significant order as 

a = \ E ((*« - m 2 lK- (O 2 ) +\ E - /?V ) 2 A S - (Q« ia2 ) 2 ) +g E o 



Q 



1 2 3 1 3 3 

77T E/ ^<y.\^o.iQo.\ai + E/ ^01^2^0102 "I" El ^ai^a^Qm a-> ~ ^7 E t aJ>a 2 Qa2a 3 Qa 1 a 3 



1 12 3 

~~ T^j y ' Q a\a 2 Q a 2 azQ a\az ^ ' t ai t a ^ aia ^ a20l ^ a30 , 4 ^| ^ ' Q a^afi aictzQ ct2Ctifd 0.30.4, 

~~ El ^<*1 Qa9.a3Qa1a9.a3 TTj" El tat tg?tanQai q^q :i El Qai 020304^^1 a9 Qazan ) 

where summations are over the unordered distinct indices. The difference in this 
expression from a simple combination of two fully or sparsely connected systems is 
in the temperature dependence of terms Ai, A2 and in the entropic terms of coefficient 
i>2 - the sparsely induced ferromagnetic/mixed solutions differs at second order. 

The simplest solution to this equation, consistent with a non-zero solution and 
the n — > limit, is a replica symmetric one plus fluctuations: taking t a = t + S a , 

Qa\a2 Q ~\~ Raia2i Q01Q203 Qz 0~Qaia2az and Qaia 2 aza4 Qa ^Qa\a2aza4 ■ These 

fluctuations should be considered independent upto permutations on the set of indices. 
We may now rewrite the auxiliary free energy (21) in terms of the RS solution and 
fluctuations as 

A = A(t, Q, Q 3 , Q 4 ) + W E S a + % E R aia 2 + ^3 E ^Qaia 2 a 3 + ^4 E $Qa 1 a 2 aza 4 (22) 

V 

+ ^E(^") 2 + ^E s a Si3+2c^R a pS a + vY J RapS 1 +-Y.Rip+QY, Ra^R ai 

+ "J E RaffR^/S + E ^ Qa0l (X5 Qa0~/ + Y R a0 + ZS a ) + E $ QaffyS (^ 4 (5Q a g 7 j + Y±R a p) . 

The prefactors are chosen to make connection with a standard stability analysis result 
which we use in the following section [dAT78]. For the RS solution to be a saddlepoint 
solution the terms linear in the fluctuations must vanish so that: 

w = o = d** + 1 ( ri + Q + \e - 2 g 2 - At 4 - Igt 2 + ^g 3 - ^(t 2 + g)) (23) 

2=0 = + \r 2 Q + g 2 - it 2 + 2gt 2 - ^g 3 + V - ^g 2 t 2 + ^ 2 g 3 t - 

z 3 = o = ^M 3 + 3g) , Z 4 = = ^ (25) 

These allow paramagnetic (t = Q = 0), ferromagnetic (t ^ 0,Q 7^ 0) and SG 
(t = 0,g 7^ 0) solutions. 



^.2. Stability near the high temperature transition 

The RS solution has t and Q such that the exponent is a local extremum. The 
saddlepoint is a stable local maxima if the quadratic form, which may be described by 
a real symmetric Hessian, is positive definite. This is difficult to calculate in the general 
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case with non-zero fluctuations in all the order parameters. Instead we wish to consider 
a restricted set of fluctuations with 5Q ai0l20l3 = 5Q aia2a3a4 = 0, thus we are testing 
coupled instabilities with respect to marginal magnetisations and spin-spin correlations 
only. This restricts somewhat the set of possible perturbations but provides a concise 
approximation to the stability properties, including the Replica Symmetry Breaking 
(RSB) ansatz, we expect to generalise very well to inclusion of more complicated 
multispin instabilities. 

We can in the restricted case calculate the eigenvalues of the quadratic form which 
are degenerate and only of three types in the limit n — > [dAT78]. The stabilities 
depend on a combination of coefficients given by 

C = C - V = t (-1 - yQ 2 + 2Q + ^-Qt 2 + || (3Q + v 2 t 2 )^) (27) 

v=r i +e + \ Q2 (28) 

Q = -|-^ + ^ 2 + 3gt 2 (29) 

K = -{^ + l )Q 2 - 2t *Q> ( 3 °) 

where we have substituted the definitions (25) for the higher order parameters. These 
combine to give the longitudinal instabilities 

\± = ±(a' + V-4Q + 311± y/(A' -V + 4Q- 3ft) 2 - 8(C) 2 ) , (31) 

and the replicon instability 

\i = V -2Q + U . (32) 

By an expansion of the coefficients, in terms of A 2 = X s — 1 near the SG boundary 
and Ai = — 1 in the vicinity of a ferromagnetic boundary, it is possible to determine 
leading order RS solutions and their stability properties. We consider the two cases that 
one component is small and the other large, and the case of a fine balance between the 
two (Ai = 0(A 2 )). Results should be interpreted with reference to figure 1. 

4-3. Solutions below the high temperature transitions 

4-3.1. Paramagnetic (P) solution. The P solution (t = 0, Q = 0) is unstable 
everywhere below the high temperature transition point and stable everywhere above 
it; eigenvalues being proportional to — Ai and — A 2 . 

4-3.2. Ferromagnetic (F) solution. In the regime where A2 < and Ai > 0, only 
the F and P solutions exist, one finds the F solution to leading order has r 2 Q = t 2 = 
3r 2 Ai/(3 + r 2 ). The longitudinal (A + ), and replicon eigenvalues are proportional to r 2 , 
which is positive in this regime. The other longitudinal eigenvalue is proportional to 
Ai, which is positive and coincident with the high temperature boundary. 
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4-3.3. Spin Glass (SG) solution. In the regime where Ai C and A 2 > only the 
SG and P solutions exist. The paramagnetic solution is unstable in a longitudinal 
mode corresponding to the high temperature boundary. One finds the SG solution gives 
Q = A 2 /2,t = to leading order. The longitudinal eigenvalue A + is coincident with 
ri, which is positive. The other longitudinal eigenvalue is proportional to A 2 , which is 
positive. The replicon instability is given to leading order as 

a -^B-^ + ^)' (33) 

which may be negative or positive depending on the value r^. The expected replica 
symmetry breaking instability is attained when J or p is large, but where p is near the 
percolation threshold and with J small the spin glass can be stable. Since our analysis 
includes the VB model as a special case, which is prooven to be unstable in the RS 
spin glass solution [Mot 8 7], this result must indicate some pathology in the stability 
analysis - the failure to consider higher order fluctuations. Nevertheless our results 
may be indicative of a general weakening of instabilities as the sparse coupling limit is 
approached. 

4-3.4- Co-emergence regime X s ~ \ F . In this case we make an expansion with both 
Ai and A 2 small, and consider both the SG and F solution. The SG solution and 
eigenvalues are unchanged at leading order except in 



(34) 



x_ = | (3, - 2 - vTTw) + a? + ^ - \ + 0(2 - „) 

where A 2 = pA x . This coefficient is negative provided < p < 2, so a longitudinal 
instability emerges when 

[A c = Ai - A 2 /2] = . (35) 

The ferromagnetic solution is changed in both the RS order parameters 

.2_oA2 AA . ^ _ A , A2 + ^ , 3W 2 \ 



t 2 = 2A -A X A 2 ; Q = A 1 + A — ^ + ^ (36) 

V 3 r 3 / 

at leading orders, and in eigenvalues 

A- = + A? g?) + A?0(2 - (37) 

^A^ + ^-AA^^) ,33, 

The other eigenvalue remains coincidence with A 1 for all p (indicating stability). Thus at 
leading order we have a transition from an RS ferromagnetic system to a longitudinally 
stable SG when p = 2 (35). 

Consider first the situation in which at the triple point r 3 is not sufficiently small 
in comparison with r 4 to make (38) positive, or that the alignment of the ferromagnetic 
moment is in the dense part t> 2 = 0. The second order corrections to the eigenvalues for 
the ferromagnetic regime indicate that in the vicinity of the transition the replicon mode 
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(38) is more negative than the non-negative valued longitudinal mode (37). We have a 
negative value for the replicon instability, but positive value for the longitudinal stability 
on the critical line, so that at second order we can predict the existence of a RSB unstable 
phase of non-zero magnetic moment (a mixed phase, M) between the ferromagnetic and 
SG solutions in a region about the critical line. This result is qualitatively similar to 
the Viana Bray (sparse) model stability result for variation of p [VB85] . If the replicon 
eigenvalue is positive on the critical line the mixed phase may be shifted marginally 
about the critical line. Unfortunately the term r 3 is sufficiently small to cause positivity 
in the second order term for many ensembles. A similar complicated dependence of 
the longitudinal mode exists for the spin glass solution 34 - the results in combination 
suggest existence of some systems with longitudinally stable F-SG coexistence regimes, 
rather than a mixed state, near the critical line. However, we suspect these details to 
be artefacts of the restricted stability analysis. 

4-4- Concluding remarks on high temperature solutions 

The RS ferromagnetic solution is the unique stable RS solution where it exists, except 
near the SG transition point (35). The SG RS solution is we suspect an unstable one in 
the replicon mode, although the stability analysis indicates that at higher temperature 
there may be a region in which it is an RS stable solution. The paramagnetic solution 
is unstable everywhere below the high temperature transition lines. The ferromagnetic 
phase becomes unstable to replica symmetry breaking (in {t a , Q ai a 2 }) in the vicinity of 
the transition, generating a mixed phase. 

It is of course essential to consider the line (35) represented in terms of some 
variation of our parameters. This line may be calculated to leading order as a function 
of £7 and 5f3, which is undertaken in Appendix C. In so doing we find the line is typically 
orientated towards 7 greater than 7c as temperature is lowered, indicating a transition 
from 'dense' to 'sparse'-type order; a result which is applicable for either type of order 
(SG or F) in the dense part. This result also implies, interestingly, that ergodicity 
breaking may disappear as temperature is lowered; a sufficient dense-coupling tendency 
towards ferromagnetism may allow such a transition in the case of weakly ordered sparse 
SG. One can consider higher order corrections to the eigenvalues in the vicinity of the 
triple-point, and in so doing we expect that in most cases a mixed solution will be found, 
that in some region of parameter space the RS ferromagnetic solution will be unstable 
when second order corrections are included. 

This analysis leaves open the corrections to the stability analysis due to higher 
order parameters, and results at second order: (34), (38) and (34), motivate such a 
consideration. A sensible way to probe higher order parameters might be to consider 
other restricted sets of fluctuations, SQ aia2a3 a function of St a and SQ ai0l2 for example. 
We have found that complicated terms in r 3 ,r 4 appear also in these cases, but we 
anticipate increased instability if fluctuations are allowed in all 4 order parameter types. 
It is likely such an analytic effort would however be better directed to a more general 
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Figure 1. A qualitative description of RS solutions and their instabilities, the auxiliary 
system results being valid in the vicinity of the high temperature transition line. Dot- 
dashed lines indicate emergence of high temperature phase transition. The darker 
(upper) lines indicate the relevant high temperature solution from this set. Left - 
The high temperature transition can be of ferromagnetic or SG type. If the types of 
order in the sparse and dense part are not complementary there can exist a triple- 
point (7c). Below the triple-point the type of order will in most cases be that 
induced by the dense subsystem. The additional dot-dashed line descending from 
the triple-point corresponds to the RS F-SG transition line (35), the RS ferromagnet 
being the stable solution between this line and the corresponding Ferromagnetic high 
temperature transition. The triplepoint analysis is valid where A^' < 1. One 
can reverse the labelling F,SG and corresponding A symbols in the diagram, with 
the qualitative nature of the diagram being similar. Centre - In cases where the 
dense distribution has negative mean coupling the ferromagnetic solution to (19) may 
disappear discontinuously, so that the high temperature transition is a discontinuous 
one, the triple-point analysis applies to no part of the diagram. Right - In the case of 
competing alignments of ferromagnetic order then there will exist competition between 
these alignments; again at a triple-point, the solution below this point can have (q a , q a ) 
non-zero in a combination determined by (14), the free-energy is symmetric about a 
dotted line (3 Jo — pT\ at leading order in A. This is a vertical line in the (3 — 7 plane 
at leading order. 



1RSB formulation solved by numerical methods. The case of a degenerate transition 
(A^ ~ A+) by a comparable method is considered in Appendix B.l. 

5. Low temperature results by RS assumption 

5.1. Replica symmetric ansatz 

The RS ansatz is here applied to the full free energy expression, the order parameters 
being invariant with respect to relabelling of spins. Returning to our order parameter 
description in terms of 7r(<x) we have the standard definitions 

tt(o-) = J dW(h) 2cogh /^ ! Qa = Qi ; q(a ia2 ) = q2 (39) 

so that W(h) is now a normalised distribution on the real line describing the moments of 
7r(er), which must be determined. The self consistent equations may then be reexpressed 
in terms of the set of parameters W, q~i, q2 at the saddlepoint 

\{h-H)\ (40) 

' 6,L,Ai,A 2 



W(h) oclf[ J dW(hi)d<f>(xi 

\i=l 
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1 L 



if = ftj,,^ + + h s + AiV + A 2 x + - XI atanh ( tanh (^) tanh(/3/ij)) (43) 

P i=i 

The average in 6 is with respect to 6, A averages are over Gaussian distributions of zero 
mean and unit variance, and L average is over a Poissonian probability distribution of 
parameter p. For 6 = one may reduce the expression (41) to a simple one, as a moment 
of the distribution W(h). 

This equation may be solved by the standard method of population dynam- 
ics [GBM01], the distribution W(h) is represented at time t by a histogram of N fields 
{hi}, alongside the scalar order parameters and . During any update step a sin- 
gle field (i) is randomly selected and updated according to a sample from the quenched 
parameters {L, b, x, h, Ai, A 2 }: 

hf +1) = b^M, + b^hP + h S + \?Jjq~ 2 + \f X 

+ W atanh(tanh(/34) tanh(/3/i* )) (44) 

P j=l 

u(t+l) _ At) . At+l) _ At) 

all other fields being left invariant. Following this we must update the other order 
parameters, the field is removed and the new field hf +V> is added to obtain 

g? +1 > = qf ) + tanh(/3/i? +1) ) - tanh^/if) 

= q t 1] + tanh 2 (/3/4 <+1) ) - tanh 2 ^) . (46) 

Note that in the case that b ^ 1 it is necessary to associate with each field in the 
histogram the sampled disorder in the relevant field update step, in order that qf 
may be incrementally updated. 

Convergence through this procedure to the correct solution, to within numerical 
accuracy dependent on N, is fairly robust. In order to avoid systematic errors, and 
attain convergence in a suitable time scale we must carefully choose initial conditions 
and population size, decide upon convergence criteria, and a sufficient level of sampling, 
as discussed in Appendix D. 

The non- variational free energy (3 (/), coincident with A after the appropriate limits 
have been taken, may be written as 

—A = - ^ (log cosh (3x) + (log(l + tanh (3x tanh /3hi tanh (3h 2 )) - (log 2 cosh(/3/i)) 

1 1 1 

+ p (log cosh atanh(tanh(/3/i) tanh(/5x))) + -pJ qf - ^P 2 J(1 ~ ~ ^ 

With appropriate scaling of the coupling distribution <fi and 7 with N it is possible to 
show equivalence of the sparse part with the dense part at large connectivity p. 
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Other quantities of interest are 

(S ai ...S ak ) = J dW(h) tanh(/3/i) , (48) 

where W(h), the auxiliary field distribution, is identical to the local field distribution. 
A sufficient statistic to describe the field distribution along the dense alignment is q~i 
combined with q 2 in the case 6 = 0. The order parameter q x is related to the correlation 
of spins, precisely as the mean spin along the alignment, 



9i = (^E 6 ^); ( 49 ) 




this definition is achieved by the derivative of the variational free energy with respect 
to hp. Similarly, the variance in the external fields, x 2 ; is conjugate to the parameter 
q 2 in the replica formalism, and h s to the mean alignment along sparse ferromagnetic 
orientation. q 2 (which may also be determined from the local field distribution (48)) 
is the SG order parameter, q~i and qi represent ferromagnetic type order parameters. 
More generally it is possible to show that the local field distribution amongst replica 
spins (48) is analogous to the field distribution for real spins in a typical sample from 
the ensemble X. 

The entropy is an additional physical quantity of interest - it is known that this can 
become negative at low temperatures in cases where the RS ansatz is insufficient and 
provides an indication for an over-simplistic assumption. The energy may be determined 
simply from the free energy as 

e= -p/2 J d<i){x)xt&nh(3x-J q 2 1 /2-(3{jq 2 2 -l)/2-h D q 1 -h s J dW{h)t&nh{(3h) 

i~ f , TTr /, x , TT ./, x ,,/ x (1 _ tanh 2 (8x)) tanh(/3/ii) tanh(/3/i 2 ) „ , ,s 
- p/2 J iW {hl )iW(W( X ) X ( - l - ^Jfe^ -fafe-l) (50) 

and the entropy calculated by the Helmholtz relation 

8 = P{e -(f)). (51) 

It is possible to see that these thermodynamic extensive quantities (51), (50) and (47) 
differ from the summation of two independent VB and SK subsystems only in the 
saddlepoint value of W(h). 



5.2. Longitudinal stability analysis 

Recent work has shown that it is also possible to test the so called longitudinal 
instabilities within the framework of population dynamics. Consider that a solution 
to the equations (41), (42) and (40) can be found, one can then examine whether the 
trajectories of fields in the population {hi} are stable against small perturbations. An 
unstable trajectory is indicative of ergodicity breaking and may be characterised by 
divergences in certain properties. It has been shown by Kabashima [Kab03] that the 
method we shall outline directly tests the longitudinal stability eigenvalue in the AT 
formalism [dAT78] for the simpler case restricted to dense couplings. In a sparse coupled 
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system Rivoire et al [RBMM04] related the divergence in the mean square fluctuation 
to the determination of SG susceptibility through the fluctuation dissipation relation. 

A more general way to test local instabilities in the sparse case is to consider a 
less restricted ansatz (1-step of replica symmetry breaking (1-RSB) [MPV87]) on the 
order parameters, and determine if the state found is identical to the simpler ansatz. 
To our knowledge there is no formalism sufficient to test all longitudinal stabilities for 
general topologies. We would expect our system to conform to the hierarchy of replica- 
symmetry broken ansatze, but possibly not be 1-RSB anywhere. We do not attempt to 
extend the current analysis to the 1-RSB formalism due to the computational difficulties, 
especially at low temperatures; moreover, we find the RS-based analysis to provide a 
good description in much of the phase space, and trust it to be indicative of trends even 
in the replica-symmetry broken phases for the simple statistics studied. 

Instability within the population dynamics solution arises out of the microscopic 
processes occurring in the update equations. In order to correctly define the 
consequences of microscopic variation it is necessary to reformulate the dense part of the 
field update (44) and represent this by the microscopic processes of iteration of a full set 
of fields on a tree (Bethe Lattice), alongside the sparse process. This is equivalent to a 
reformulation of the problem as belief propagation [Kab03], and testing the convergence 
of the dynamics to a unique solution. 

It is necessary to represent the dense interactions in the saddlepoint equation (46) 
by a similar structure to the sparse interactions, as a summation over many spin-spin 
interactions rather than their statistical average. Here, the saddle point equations of 
(46) take a microscopic form 

bpq x + \fq 2 -> i &S* +1) atanh (tanh(/3^' ) ) tanh hf) , (52) 

' j 

where are now random samples from the dense coupling distribution (which may 
be taken as a Gaussian jV(Jq/N, J/N)), and bf +1 ^ is the orientation of the spin on 
the target site, j runs over all fields excluding those involved in the sparse subsystem. 
We will now calculate how a small variation in the incoming fields acts to perturb the 
outgoing field in a particular update. 

hf +1) + = J2 ^-atanh (tanh(/3af ) tanh(/if + 5hf)) 

j\{i,il-i L } r> 

+ \Y. atanh (tanh^af) tanh^f + Shg)) (53) 
P j= i 

Where L, the set of fields {ij}, the sparse and dense couplings and the target disorder 
are random samples in any update. Assuming the perturbations 5h to be small and 
knowing the couplings x we can make a Taylor expansion to linear order obtaining 

e" e (i-t^w^ay+i: V < Tt ( 5ll!"tS ) < - (54) 

j\{i,h...i L } j=i 1 - tanh(/3x} ; ) tanh(/ij/) 
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with the indices again referring to the same quenched samples. Using the fact that the 
first part is a sum of N — as iV random variables we can simplify the description of 
this term by the central limit theorem, describing the processes by a normal distribution 
of mean and variance 



which characterise linear and non-linear (SG) susceptibilities by analogy with results 
for sparse [RBMM04] and dense [Kab03] systems. We expect a negative value in xi 
to indicate local linear stability, while a negative value in X2 indicates local non-linear 
stability. We expect the non-linear stability to be sufficient to detect ergodicity breaking 
trends in solutions. 

In the case of a paramagnetic solution results may be determined exactly, the 
linear instability in 5h is found to correspond to the ferromagnetic solution criteria 
(19). The non-linear susceptibility is found to coincide with the SG solution criteria 
(16). Therefore the paramagnetic solution instabilities determined in section 4.2 are 
reproduced exactly by this method. 

Generally we note that by the nature of the population dynamics algorithm there 
exist a number of noisy effects: (i) Numerical finite precision errors; (ii) The noise arising 
from the random order in which fields are updated; (iii) the systematic effect of updating 
first the field hf \ then qf 1 and q®; (iv) the finite precision in the histogram; and most 
importantly (v) the random sampling of quenched parameters in each field update. We 
can therefore expect convergence of our dynamics to limit not only the final precision 
achieved, but also the set of solutions, to those stable against such perturbations. We 
assume this class of perturbations to be sufficient to probe the linear stability even in the 
low temperature regime. Any linear instability should be observable in some moment of 
the histogram so provided we test convergence in W effectively we hope to rule 

out linear instabilities. 

For the other phases the instabilities depend on the details of the field distribution. 
Consider that we have arrived at some fixed point described by {^}- We can test the 

stability of this point by adding small perturbations {(5/if^) 2 } and determining \2 once 
the dominate modes emerge. Details of the stability analysis implentation are given 
in Appendix D. 

6. Numerical results of population dynamics 



/3J ((1 -tfwh/3hi)6hi)/N ; 1 J ((1 - tanh/^) 2 ^) /N , (55) 



respectively. 

We are interested in two types of quantities for any solution 





(56) 



We employed population dynamics to investigate a number of systems concentrating on 
those in which we could examine the effects of competing order as both a weak effect 
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(near the high temperature boundary) and a strong effect at lower temperatures. Many 
critical quantities could not be sufficiently resolved at high temperature primarily for 
reasons of finite precision, and near the percolation threshold, in our algorithms we 
present cases with (f3//3 c < 10) and p = 2, 1//3 C being the largest critical temperature 
(typically 1 for examples chosen). 

6.1. Data format 

All data is based on linear spline interpolation of array data. In the figures over the range 
7 = (0, 1) we use a point spacing of (5/3, 5j) = (0.025,0.025). For these diagrams we 
present results based on a single run from random initial conditions, and plot the mean 
and error bars (linear spline interpolation based on values at each sample point) over 
20 samples of data. These 20 samples are selected in successive time steps (a time step 
is iV single field updates, 1 for each field), following the convergence of the distribution. 
Therefore, the error bars are not over independent samples, and as such not necessarily 
well described by an uncorrelated Gaussian distribution, but this assumption gives a 
first approximation to single time-step fluctuations in the measured state. Each point 
that did not converge within 500 time steps is marked by a cross. 

Figures that focus on part of the 7 range use point spacing of (0.005,0.005). We 
average quantities in these figures over 10 runs based on different seeds for the random 
number generator (different initial conditions and updates). The data point is taken to 
be the mean of 20 samples from each run. Unlike the sampling within a single run (which 
may be correlated in successive update steps) we can be confident of the independence of 
these samples conditional on the parameterisation and convergence criteria, and present 
the mean and error bars for quantities of interest (again interpolation). 

Population dynamics is not required in the paramagnetic region, for which 
simple boundary conditions based on exact knowledge are imposed, this creates some 
unevenness in certain quantities very close to the boundary. A second source of 
unevenness is the array like nature of the data points. Finally our convergence criteria 
appears not strict enough to prevent some systematic errors (drift) near critical regions, 
otherwise we expect error bars to be representative. 

Cases of incomplete convergence are treated for purposes of data collection and 
interpolation as converged results. Incomplete convergence occurs: (i) in a small number 
of cases close to critical transitions, (ii) in the case of competing ferromagnetic alignment. 
After the maximum number of population iterations is completed (500) we select the 
histogram (amongst 4 differing in initial conditions) of minimum free energy which is 
good in case (ii), but does not resolve case (i). However, we found experimentally that in 
cases (i) many quantities did not vary greatly in absolute terms between the histograms 
of different initial conditions, except for some systematic drift in the stability parameter 
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6.2. Figures presented and coupling distributions considered 

As there is a large number of composite systems that could be considered and analysed, 
it would be useful to review the choice of the specific cases presented here. 

We present data for the several aligned (6=1) combinations of sparse and dense 
subsystems. A mixture of a dense subsystem of ferromagnetic couplings with a sparse 
SG (figure 2) or anti-ferromagnetic (figure 4) couplings; and the converse, a subsystem 
of sparse ferromagnetic couplings mixed with a dense subsystem of SG (figure 3) or anti- 
ferromagnetic (figures 5 and 6) couplings. We also consider the case of two subsystems 
both with ferromagnetic couplings but unaligned (b = 0, figures 7 and 8). 

The dense couplings are parameterised by (Jo, J), so that for example a dense anti- 
ferromagnetic system is (—1, 0) in this notation. The sparse system we choose is the ±J 
model, a canonical case convenient for numerical reasons. In this model is bimodal 
and described by p the probability of sampling +J as opposed to — J for each bond. 
Therefore the sparse part is parameterised by (p, J) with p = 2, describing the relevant 
subsystems we present. 

Owing to the large variety of parameterisations we are not able to test sufficient 
systems to be able to make so general statements as are possible in the exact analysis 
of the high temperature behaviour. However, we believe these systems to be the most 
intuitive of combinations. We examined the effect of changing the distribution of sparse 
couplings (j) from ±J to Gaussian but found only marginal variations. 

6.3. Stability Results 

Results are generally characterised by a dominance of the dense system below a critical 
mixing parameter 7 C which is replaced by a dominance of the sparse system at larger 
7, although the critical value and the specific properties depend on the system studied 
and the temperature value. 

We find that everywhere close to the boundary the longitudinal stability eigenvalue 
is positive, thus the RS solutions appear to be locally stable, including in the SG phase. 
The small gap observed for the SG solution between the numerical results obtained 
from population dynamics and the high temperature transition line is we suspect due to 
finite size effects in a combination of population size and number of updates. Otherwise 
the SG solutions are unstable to ergodicity breaking as expected. The RS ferromagnetic 
solution may also become unstable at sufficiently low temperatures. Therefore it appears 
that there is a mixed phase as predicted by the stability analysis in which the magnetic 
moment is non-zero but ergodicity breaking is present. Dense ferromagnetic phases 
appear less susceptible to this instability. 
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6.4- The Co- emergence regime: properties on lowering temperature about a P-F-SG 
triple-point, 7 ps 7^ 

We can examine the type of order observed as we decrease the temperature from the 
triple-point, which provides information as to the relative importance of couplings in 
weakly ordered systems. It appears that in the vicinity of the triple point one can lower 
the temperature and finds that the type of order present is closer to that induced by 
the dense couplings in most cases (inline with the predictions of the high temperature 
analysis). The sparse type order parameter in some cases completely disappears as 
temperature is lowered. 

If the dense type order is ferromagnetic (figure 2) it appears that by lowering 
temperature a region of ergodicity breaking may be encountered, but the magnetic 
moment does not disappear; such a scenario describes a mixed phase. Conversely if one 
lowers the temperature where the high temperature is a sparsely induced ferromagnetic 
one (figure 3), it is possible not only for the ferromagnet to become unstable towards 
a mixed phase, but for the magnetic moment to disappear entirely. This characterises 
a standard reentrant behaviour, as the magnetic moment re-enters the value it had in 
the paramagnetic phase. For the anti-ferromagnetic SG state (figure 5) there are some 
indications of transitions from SG to mixed as well as ferromagnetic to mixed transitions 
as temperature decreases. 

A more exotic possibility observed in the simple combination of a dense ferromagnet 
and sparse SG subsystems (figures 2 and 4) is that a SG phase may become an RS 
ferromagnet and then a mixed phase as temperature is lowered. The numerical results 
are inconclusive on this point, but one can combine this with knowledge of the exact 
transition line result at the boundary (35), and the intuition that ergodicity breaking is 
more likely at low temperatures in frustrated systems. It is unusual to our knowledge for 
ergodicity breaking to disappear as temperature is lowered, only for it then to reappear. 

6.5. Discontinuous high temperature transitions 

Figures 5 and 6 consider scenarios in which the high temperature transition is potentially 
discontinuous. Where the boundary is discontinuous, the sparse coupling order is 
ferromagnetic; it appears that the sparse-induced ferromagnetic order overwhelms any 
emergence of SG order in the vicinity of the discontinuity (7 7c), so that the SG phase 
must occupy only some very narrow region about the high temperature transition. We 
can anticipate for both the SG and paramagnetic phases, that convergence of population 
dynamics might be subject to especially strong finite size effects in the vicinity of this 
transition. We vary the sparse ensemble through J to generate the different scenarios, 
p may also be varied to create this effect. 

With the transition nearly discontinuous (as is the case for J = 1/2 figure 5 in which 
the discontinuity is visually imperceptable) one observes a similar behaviour to the case 
of a SG-F combination (figure 2) except in the weaker suppression of the ferromagnetic 
order at low temperatures. Clearly, in this case the antiferromagnetic dense component 
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cannot induce any order except for a paramagnetic one and the emergence of a SG is a 
property of the sparse distribution. 

In figure 6a where a discontinuity is clearly visible we find ergodicity breaking is 
surpressed in the vicinity of the discontinuity. The ferromagnetic phase dominates near 
7 C at high temperature but gives way to a mixed phase at lower temperature. Figure 
6b shows that in the continuous case the high temperature prediction is qualitatively 
accurate. The SG state appears to dominate near j c , with mixed phases appearing at 
lower temperature. 




Figure 2. A dense subsystem (Jo = 1, J = 0) with a sparse subsystem (p = 0.5, J = 1) 
- dense ferromagnetic and sparse SG couplings. A schematic denotes our qualitative 
understanding of the solutions. Main figure: The exact critical lines ((16), (17), (19)) 
are indicated by dot dashed lines with the darker (upper) line being the exact high 
temperature transition. The thick line \2 — 0, and the line of zero magnetisation 
qi = 0, indicate transition points as determined by the population dynamics algorithm. 
The continuation of q\ =0 must be inferred by the contours in q\ (thin lines). 
Error bars in all contour lines are plotted and indicate time-step fluctuations within a 
single run. Transitions are between unstable-stable and SG-F solutions, the unstable 
ferromagnetic solution is classified a mixed state. It appears the effect of sparse 
frustrating bonds is to make the RS dense ferromagnet unstable as f3 increases. 
Similarly the sparsely induced SG may be susceptible to a transition towards a mixed 
state with increasing (3. The mixed phase remains close to 7 = j c , thus the competition 
of interaction types remains in some sense balanced as temperature decreases. Inset 
figure: The right figure shows the same contours as the main figure in magnification 
with corresponding error bars (variation between different independent minimisations). 
We also plot the high temperature critical line (35), which indicates the exact gradient 
for the curves xi = and q\ — 0, this appears not to be closely followed by the line 
Xi which is due to finite size numerical errors and coarseness of point sampling. 
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Figure 3. A dense subsystem (Jo = 0, J = 1) with a sparse subsystem (p = 1, J = i), 
i.e., dense SG and sparse uniform ferromagnetic couplings. A schematic denotes our 
qualitative understanding of the solutions. Main figure: Lines and symbols are the 
same as for figure 2, in addition the thick dotted line represents a lower temperature 
bound to positive entropy (include the negligible time-step errors), below this line 
the method is inconsistent as would be suggested by our results on instability of the 
SG phase (to which this curve belongs). Inset: The low magnetisation contours, and 
stability parameter follow quite closely the triple-point analysis prediction (35) in the 
vicinity of the boundary (inset). These lines break indicating a mixed phase. It appears 
there may be no transitions in the SG state towards, a mixed phase. The mixed phase 
is shifted to higher 7 with increasing /3 indicating the greater relative importance 
of the dense disordered couplings in determining phase behaviour as temperature is 
lowered. The ferromagnetic state is first unstable to a mixed state, which may in turn 
be unstable towards a SG state. 

6.6. Unaligned ferromagnetic couplings 6 = 

This is an interesting case of competing ferromagnetic alignments; as one decreases the 
temperature in the vicinity of 7c it appears that, as a thermodynamic solution, the 
sparse alignment is marginally dominated in the thermodynamic state by the dense 
alignment (figure 7). However, the sparse aligned solution persists at lower temperature 
as a metastable state, responsible for the non- convergence, which is locally stable against 
ergodicity breaking (% 2 < 0). 

In fact, metastable states of both dense and sparse alignment exist across a 
wide range of parameters 7 for any (3 below the triple point (figure 8), indicated by 
populations converging to different numerically stable solutions. These metastable states 
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Figure 4. The combination of dense subsystem (Jo = 1, J = 0) and sparse subsystem 
(p = 0, J = j) is a combination of ferromagnetic and antiferromagnetic (SG) effects; 
lines and symbols are as in the previous figures. Not surprisingly this case is almost 
identical to the case of figure 2. There is excellent agreement between the thin data lines 
(ferromagnetic moments) and the high temperature prediction for their disappearance 
(inset). A clear cusp exists in the stability, indicating that an ergodic ferromagnetic 
behaviour may emerge in some systems as temperature is lowered, this indicates 
increasing importance of the dense couplings as temperature is lowered. However, 
instabilities do reemerge at lower temperature. The SG is susceptible to a mixed state 
as temperature is lowered. 

are characterised numerically by q~i > 0, qi ~ and qi > 0, q~i ~ and we speculate 
the thermodynamic state involves a first order transition between the two solutions. A 
necessary condition for the existence of a metastable state is the non-negativity of both 
A± — 1, but this is not necessarily sufficient. 

Histograms with a well defined alignment tend to converge to the locally stable 
state with the same corresponding alignment, whether or not this is the thermodynamic 
solution. Random initial conditions appear to converge towards the thermodynamic 
solution at high temperature, but have a slight bias towards the dense alignmnent at 
lower temperatures (atleast over a small number of updates). It is difficult to determine 
if the solution space generate a larger basin of attraction for the dense metastable 
state, or one, in some sense, proportional to the free energy. This may be important 
given that we envisage composite systems as components in optimisation problems were 
dynamical minimisation by local search, similar to population dynamics, is a critical 
feature, possibly more so than any static properties of the ensemble. 
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Figure 5. The combination of dense subsystem (J = — 1, J = 0) and sparse 
subsystem (p = 1, J = h) is a combination of antiferromagnetic and ferromagnetic 
effects; lines and symbols as in the previous figures. The dense couplings cannot 
independently induce SG or ferromagnetic order, but do suppress ferromagnetic order. 
The resulting behaviour is comparable to a SG ferromagnet combination- figure 3. By 
contrast to the dense SG coupling there is weaker suppression of the magnetic moment 
at low temperature, and ergodicity breaking is a weaker effect at low temperature 
in the ferromagnetic phase. In this diagram the high temperature transition line are 
marginally discontinuous as discussed (section 6.5). Therefore we do not provide a 
detailed inset for the triple point analysis region. 



7. Conclusion 

This paper considers the question of systems consisting of a combination of sparse 
and dense random couplings. Given the complementary properties of these types 
of structures, combined with the possibility of a robust statistical description, many 
engineered systems may benefit by a decomposition as laid out here. Such an example 
we are currently pursuing involves the combination of sparse and dense spreading codes 
for multiuser access [RS07]. An alternative application of this analysis might be in 
modelling random attacks on network structures, where correlations between attacked 
elements are induced [HM]. 

This paper has considered a simple model in which canonical sparse and dense 
disordered models have been mixed. In so doing we have observed some clear trends 
which may generalise. These include the tendency for dense induced order to dominate 
sparse induced order when in competition. In a system with a dense ferromagnetic 
tendency, it is possible for this order to emerge from a sparse induced SG as a replica- 
symmetric stable system by lowering the temperature, provided the order is sufficiently 
weak (high temperature). In the vicinity of a high temperature transition the effect 
of suppressing ferromagnetic order by an increasing antiferromagnetic tendency in the 
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Figure 6. Anti-ferromagnetic couplings in the dense part aligned with ferromagnetic 
order in sparse couplings can cause the high temperature ferromagnetic solution criteria 
(18) to be unmet for any (3 in some range of 7. If 19 describes the high temperature 
transition as 7 decreases towards this range then there is a discontinuity (a) and the 
stability analysis for the region A 2 = 0(Ai) does not apply, otherwise the transition is 
a continuous one (b) with the stability analysis applicable, (a) With antiferromagnetic 
dense couplings (Jo = —2, J = 0) and ferromagnetic sparse couplings (jp — 1, J = 1), 
the transition is discontinuous. The sparse order dominates the dense order near 
7 w 7c, the relative importance of the dense couplings increases as temperature is 
lowered. The closeness of the magnetisation contours to the line X s = 1 may be 
a significant finite size effect, these lines do not coincide in the thermodynamic limit, 
(b) With antiferromagnetic dense couplings (Jo = — 4, J = 0) and ferromagnetic sparse 
couplings (p — 1,J — 1), the transition is continuous. Results appear to follow closely 
the high temperature prediction. At lower temperature 7 ~ 7c appears to remain 
approximately coincident with a SG-M transition, the ferromagnet is again unstable 
to a mixed phase in many systems. 

bonds may have a very similar effect to the introduction of frustrating effects in the 
interaction structure. Finally, the case of the unaligned system indicates that the 
dominance of dense ferromagnetic over sparse ferromagnetic couplings is marginal in 
equilibrium properties, however it appears a metastable behaviours may have non-trivial 
consequences for dynamics. 

Several questions remain open - importantly, the question of the low temperature 
limit is not resolved in this paper nor is the behaviour in the vicinity of the percolation 
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Figure 7. A system with 6=1, and uniform couplings of (Jo = 1,J = 0) 
and (p = g , J = 1 ) , models the competition between two unaligned ferromagnetic 
orderings; lines and symbols as in the previous figures. The converged solutions were 
found to be ergodic everywhere (\2 < 0) and were of positive entropy. In a part of 
the space below the cusp in the high temperature transition (<y c ) population dynamics 
converge to one of two differently aligned locally stable states as can be seen by the 
crosses. Varying 7 through the centre of this region there is a rapid decline in q~i with 
the magnetic moment growing instead in q\. We speculate that this is a first order 
transition between two locally stable states. 

threshold. These may be analysed by related statistical physics methods. There are a 
number of other variations which would allow the insight gained from the current study 
to be strengthened, including reformulations of the basic model so that the arbitrary 
parameter 7 may be shown to be only a technicality; so that related models, for example, 
in which a dense system is fixed and one gradually adds more bonds (a variation of p), 
may be surmised. One can consider a number of ways of scaling <J),Jq and J with 7 
but we hope the properties we have reported here are robust against the most sensible 
alternatives such as scaling all coupling moments of each subsystem uniformly with 7. 

Finally, with the increasing interest in complex systems of varying connectivities 
and interaction strengths, we believe the current study of a simple composite system 
comprising elements which have both sparse (and strong) and dense (weak) interactions, 
represents a first step in the principled analysis of their equilibrium behaviour. 
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Figure 8. A system with 6=1, and uniform couplings of (Jo = 1,J = 0) 
and (p = i , J = 1 ) , models the competition between two unaligned ferromagnetic 
orderings. These figures show magnifications of figure 7 near j c . Left figure: In 
the vicinity of j c a large number of states fail to converge after maximum algorithm 
time of 500 steps (dots). We then have to chose which of the 4 replica histograms 
evolved in the dynamics produces a result of smallest free energy. If we evaluate the 
magnetisation in this state we find a transition indicated by change of sign in q~\ — q± , 
this is plotted with error bars (time step fluctuations) for a single run. Left figure: 
Sampling 20 systems with a maximum time of 100 time steps in each iteration we show 
the contour q~\ — q\ = with sampling errors for several different histograms evolved 
in parallel: a sparsely aligned ferromagnetic (left most line), random SG or nearly 
paramagnetic (central thin lines) or densely aligned ferromagnetic initial histogram. 
The central thicker black line is the equilibrium result as determined by from the 
histogram of minimum free energy. With only 100 iterations there is a slightly broader 
set of unconverged states (dots), but we already see the emergence of two well defined 
dynamical transitions - left and right lines, matching approximately the result in the 
left figure. The alignments of the random histograms appear to approximately match 
the thcrmodynamical transition at high temperature but are noisy and slightly biased 
towards a dense alignment at lower temperature. 
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Appendix A. Replica calculation 

A particular ensemble X is described by a set of order one scalars j Jo, J, p, hs, hp, X 2 } 
and a sparse coupling distribution <fi of finite moments on the real line, without finite 
measure at zero. These parameters contain an unstated dependence on the mixing 
parameter, we choose this to be in the couplings (f),J and J but leave this unstated for 
brevity. An instance of the ensemble is constructed according to the following probability 
distributions, subscripts being the relevant random variables 
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The dense coupling distribution <\P may be taken as a Gaussian of mean Jq/N and 
variance J/N. The field distribution for each site is a Gaussian of mean hoh + hs and 
variance x 2 - 

We describe the properties of a particular ensemble through self averaged quantities 
calculated from the mean free energy density (self averaging assumption) which with 
the replica identity [MPV87] may be written 

(/> = lim — lim (Z n ) , 

with (5 the inverse temperature. The replicated partition function, averaged over samples 
is given by 

S 1 ...S n \ \ I a <«) J / J D \ I Q fe> J / A,J S 

* M?s>*}) h ) 6 ■ 

As is standard in such calculations [MPV87] we maintain only those terms in the 
exponent which are order N, taking asymptotic values in n — > for brevity in all 
prefactors. The quenched average on dense couplings may be factorised and taken by 
linearising the exponent, resulting in 



On composite systems 29 
The sparse coupling average may also be taken by standard sparse methods [Mon98] 
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and the field average 
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Introducing the following order parameters, $: 

7r ( <T ) = ^ 2 <Wo" ^ = TV S 6 ^ ' 9<ai"2> = ^ E ^i^ 2 ' 

j i i 

we are able to write the replicated partition function as an integral in the set of order 
parameters 

(Z n ) oc [d$ E I($)exp{-iVG' 1 } ; 

(T 1 ...(T N 
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where the function I is an indicator function for the order parameter definitions 
Appendix A. We note here that there is potentially some redundancy in the definition 
of order parameters, this is allowed for a concise and general expression. 

Representing each definition in the indicator function by a Fourier transform, 
introducing conjugate integration variables denoted with a hat, 

E I($)= fd^exp{-N{G 2 + G 3 )} ; 

(T 1 ...(T N 

G2 = -E^fa- E ?(«i«2)9(«i«2> - E 7r ( cr )^( cr ) ; 

a (aia 2 ) f 

G 3 = - log]T (exp J -&£ L° a - E ?<«i« 2 >* ai <7 aa " 

a \ { a («i« 2 ) ) 1 b 

assuming N is the suitable scaling for the conjugate variables. We have finally a site 
factorised saddlepoint form (5). 

Appendix B. Calculation details for auxiliary system 

In order to construct the auxiliary system we must determine the significance of 
terms. We keep terms sufficient to allow a second order description in A, which is 
the perturbation away from A5 = 1 and/or A+ = 1, and to allow a second order stability 
analysis, restricted to the region where A^ < 1 and A^ < 1. We calculate the expansion 
of the complicated entropic term log Tr exp {X} in equation (11) in the reduced set of 
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order parameters {t a , Q aia2 , Q ai a 2 a 3 , Qa 1 a 2 aiai\ an d to an order fit for purpose. This 
becomes for b — 1 (v 2 = 1) or b = (v 2 = 1 or 0) 

-logTrexpA" = - ^ log cosh t a -^ log cosh Q{ aia2 ) (B.l) 

- l0 § COsh Q<aia 2 a 3 > - E l0 § COsh <2<aia 2 a 3 a 4 > 
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1 + tanh(Q (aia2) ) Jl^i 

i=l 
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n 
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the sums and products are ordered in the various terms and are over the corresponding 

replica indices. The case of an aligned system and unaligned have simple averages, 

being distinguished by v 2 - The Trace requires a graphical expansion, otherwise various 

non-linear terms must in some cases be expanded upto third order. The expansion of 

(B.l) then gives 
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where sums are now over all sets of indices without repeated indices. Finally introducing 
the energetic part we have (21). 



Appendix B.l. Stability analysis for unaligned systems 6 = 

We have proposed that, given a sufficient gap between A^ and A^, then one alignment 
of the ferromagnetic solution may be considered dominant and the absolute value 
and fluctuations in the converse direction may be ignored. However, if these values 
are comparable one must consider an expansion in an additional order parameters to 
understand the ferromagnetic phase: t a — > {q a ,Qa}- To correct the auxiliary system 
we make the substitution of the type t a l = ((q a + bq a ) % )b for i th order terms in the 
entropic part (B.2). There are only two additional terms at A th and 5 th order in the free 
energy to consider, which couple the two parameters. This is a complicated expression 
to evaluate. 

Under RS we can anticipate a set of results including {q a 7^ 0, q a = 0} and 
{q a = 0, q a 7^ 0} to dominate based on numerical findings and which are both locally 
stable (section 6.6). One simple observation is that in the expansion we find that at 
leading orders, excluding critical points, there is symmetry about the line (3Jq = pT\. 
By expression of this equality in the 7, (5 plane (Appendix C) we can observe this is 
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skewed (depending on the scaling, e.g. only at subleading order (3)) towards higher 7 
with decreasing temperature. If the alignment changes monotonically with 7 then the 
point immediately below 7 C will be unbiased in alignment at leading order, but favour a 
dense alignment when the non-linear dependence is included. This qualitative argument 
appears consistent with results, but the discontinuous nature of the transition requires 
a consideration of more complicated effects. 



Appendix C. Triple point analysis 

To discriminate between the effects of different subsystems one must make an expansion 
of Aj as a function of 5(3 and £7 about the triple-point. If we take A7 = /iAT we can 
attain an equation for the critical line (35) dependent on the ensemble details 

Ai = ~^{^(3) + ^5(3 + 0(5(3 2 ) (C.l) 
f)\ s c)\ s 

A 2 = ~^{^(3) + ^5(3 + 0(5(3 2 ) . (C.2) 
Which gives an equation for the value in which ferromagnetic instability emerges 



/3J -/3 2 J+p(/3x(l-tanh(/3x)+tanh 2 (/3a;)-tanh 3 (/3a;))} 



-/3J '+ / 3 2 J72+p(/3x(l-tanh(/3x)+tanh 2 (/3a;)-tanh 3 (/3a;))} , 



4,'(x) 



(C.3) 

Pcic 



The derivatives in the denominator in J , J, and distribution 0, denoted by ' are with 
respect to 7. These are taken straightforwardly with linear scaling of couplins (3). The 
expression is quite complicated to interpret, but may be simplified further using the 
criteria for triple point existence (19) and (16). These expression may be used to find 
a triple point for a particular system in variation of some parameters (7, (3) and we can 
then test whether the dense or spin induced order dominates as temperature is lowered. 

Rewriting the transition line for linear scaling of the Hamiltonian (3), and 
resubsituting (19) and (16) we find for the aligned system (b — 1): 

(7-7c)=M/?-/fc) (C4) 
1 ~(D) + (S) 
f3 (D) /(l - 7) + (S) h 

D = tanh(/3 7 :r) - tanh 2 ((3-fx) (C.6) 

S = (3-fx(l - tanh^a;) + tanh^rr) - tanh 3 (/3 7 a;)) (C.7) 

Where the distribution (f>(x), over which averages occur, is by the decomposition (3) 
now a distribution independent of 7. Interestingly, this value is independent of the 
connectivity p, thus percolation properties appear not to be an important effect in the 
RS order stability in the vicinity of the triple-point. This connectivity independence is 
true for more general scaling scenarios than (3); whenever the derivatives of Jo, J with 
respect to 7 are proportional to Jo, J, one can make substitutions to write the expression 
as functions of only (f), 7, /3. For b ^ 1, this is not necessarily possible though certain 
properties are similar. The expressions S and D are shown in figure (CI), it is clear that 



(C.5) 
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if is negative in mean or has large variance then for most distributions (D) and (S) 
are negative, with (D) < (S). If the distribution is positive in mean and of relatively 
small variance then (D) and (S) will be positive, with (D) the larger. However, this is 
not a general result and exceptions may be constructed. 




Figure CI. The dashed lines indicate the quantity D (C.6) and the solid lines the quantity £ (C.7). 
If one has a particular distribution of sparse couplings 4> then at the triple-point one can integrate 
over these distributions to determine the quantities and hence the nature of reentrant behaviour - 
whether the weakly ordered systems are dominated by dense or sparse type order, we find the former 
in canonical cases. If the variance of the distributions is sufficient the integral for both quantities is 
negative. If the distribution is of positive mean with low variance both quantities will be positive. It is 
possible to generate many distributions for which the above generalisations do not hold, and reentrant 
behaviour of an alternative type may occur. 



Appendix D. Population dynamics numerical implementation 

The population dynamics is implemented such that the fields are assigned randomly 
at time t = in three histograms (real replica) simultaneously. We consider the 
development from initial conditions corresponding to frozen ferromagnetic, SG and 
near paramagnetic (high temperature) initial condition histograms. The corresponding 
Gaussian distributions are 

P(hf ] ) = JV(1000, 0) , JV(0, 1) and Af(0, 0.001) , (D.l) 

We also add an additional (fourth) histogram with ferromagnetic order along the dense 
alignment when 6^1, and initialise bf^ to ±1 with probabilities determined by b. 

Since the representative low-temperature order parameters are overestimated in the 
former two, and underestimated in the latter, we can be confident a solution converged 
to by all three will not be systematically biased. However, in the case of multiple locally 
stable solutions these replica remain disparate, but cover a representative set of solutions 
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quite well (e.g. 5 = 0). This parallelism is at the cost of runtime, but we found drift, or 
metastability, to be prevalent effects justifying the method. 

It is also useful to consider fluctuations of the solutions within a single histogram 
in successive time steps (N field updates). We judge convergence in both cases through 
consideration of the order parameters q±, q~i, and q 2 . This is insufficient for a system 
involving sparse interactions since the distribution is then not Gaussian, but appears 
nevertheless to be quite robust and comparable to other heuristic methods attempted. 
The criteria for termination of the algorithm is that the distributions based on all 3 (4) 
initial conditions and between successive updates converges in q±,qi and q 2 to within a 
worst case absolute precision of 0.01. Finite size fluctuations and runtime constraints 
prevent a signigicantly stronger precision. 

The field updates of (44) are random and non- sequential - but ensuring that in 
each time step, each field is updated exactly once. 

We found in almost all cases that histograms converged to a single unique 
distribution (upto finite size errors), and that this convergence was robust. When 
convergence criteria is met we analyse the histogram (replica) of lowest free energy as 
representative. Measurables are determined following convergence based on 20 samples 
over 20 time steps. For certain systems we ran upto 20 independent runs (different 
pseudo random update sequences) to check the properties were unique and fully explore 
the space. 

In the case of non- convergence by the moments criteria we halt the algorithm after 
500 time steps and consider the histogram that corresponds to the lowest free energy 
(amongst the different initial condition histograms), calculating quantities based on this 
choice. In the case that metastable solutions exists this is sufficient to identify the 
histogram corresponding to a dominant state. This occured only for certain systems 
with opposite ferromagnetic alignment (6 = 0). In all other cases it appeared solutions 
are in the vicinity of the correct saddlepoint but converging with different bises; thus 
identifying the solution of lowest free energy does not guarantee improved measurement 
of the statistics in this case. However, we judge the statistics collected to be close in 
most cases. 

Appendix D.l. Stability 

We associate with each field in the population dynamics a perturbation 5h 2 , these are 
initially chosen with a Gaussian distribution independent of the field. With 7 > 
convergence to a non-Gaussian joint distribution of perturbations and fields must be 
considered, we observed this to occur very quickly as can be determined in correlation 
functions and kurtosis, e.g. (o~(hl) 2 (hl) 2 ) c and (Shj), but could not produce a robust 
measure of convergence of the fluctuation distribution. Following convergence in the 
order parameters, we introduce the field and following 10 time steps collect data over 20 
time steps. Distribution of perturbations are carefully renormalised throughout analysis 
to account for the non-sequential updates, and the renormalisation after each time-step 
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provides an estimate for \2- Certain other numerical hacks must be brought to bear to 
prevent divergences in the cases where X2 is n °t close to 0. 
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